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Simulation  modeling  requires  accurate  input  analysis  to  ensure  validity 
of  the  study.  Hence,  the  mantra  “garbage  in  =  garbage  out.”  Much  of  the  re¬ 
search  and  simulation  code  that  has  been  written  to  date  has  been  focused  on 
traditional  parametric  methods.  Here  we  investigate  Bayesian  nonparametric 
methods  for  input  modeling  and  reliability  analysis.  Bayesian  nonparametric 
methods  have  been  shown  in  many  cases  to  produce  better  predictive  models. 
Also,  for  use  in  a  Bayesian  setting,  we  have  written  C++  classes  for  random 
variate  generation.  These  contain  functions  for  standard  and  truncated  distri¬ 
butions  as  well  as  functions  for  statistical  data  handling.  Although  we  have 
written  the  code  for  Bayesian  algorithms,  the  functions  can  be  used  anywhere 
a  good  source  of  random  variates  is  needed.  Included  is  a  detailed  description 
of  class  implementation  and  usage  along  with  complete  source  code. 
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Chapter  1 


Introduction 


Many  real-world  systems  of  interest  to  operations  research  analysts  have 
inherently  stochastic  components.  These  include,  but  are  certainly  not  limited 
to,  financial  markets,  mechanical  systems  with  uncertain  lifetimes,  and  queuing 
systems.  To  gain  a  better  understanding  of  the  underlying  mechanisms  of 
such  systems,  analysts  build  statistical  simulation  models.  The  focus  of  this 
report  is  on  input  analysis  for  statistical  simulation  modeling.  Input  analysis 
is  extremely  important  to  simulation  modeling  as  wrong  input  distributions 
will  result  in  misleading  simulation  output  and  incorrect  decisions. 

In  simulation  input  analysis  we  must  first  develop  accurate  probability 
models  to  represent  real-world  random  phenomena,  then  simulate  draws  from 
these  predictive  distributions  for  input  to  the  simulation  model.  To  date,  the 
primary  focus  for  input  analysis  has  been  on  traditional,  parametric  statistical 
methods  [see  Law  and  Kelton,  2000].  Here  we  discuss  an  alternative  method 
for  input  modeling  in  a  reliability  setting,  Bayesian  nonparametrics,  and  a 
C++  class  library  that  we’ve  developed  primarily  for  Bayesian  simulation. 

We’ve  chosen  to  investigate  Bayesian  nonparametric  methods  for  input 
modeling  in  an  attempt  to  make  data  models  more  robust  and,  as  a  result, 


1 


make  more  accurate  predictions.  With  a  Bayesian  nonparametric  model  we  do 
not  force  in  any  undue  assumptions  about  the  characteristics  of  the  underlying 
distribution,  yet  we  can  inject  prior  notions  about  the  location  and  spread  of 
the  data.  Bayesian  models  also  have  a  nice  feature  in  that  we  can  “learn  from 
experience” .  We  can  continually  update  the  model  with  new  data  and  therefore 
capture  more  and  more  information  about  the  underlying  mechanism. 

The  focus  here  is  on  a  Bayesian  nonparametric  technique  to  simulate 
cumulative  hazard  functions  using  Beta  stochastic  processes  [see  Hjort,  1990]. 
Cumulative  hazard  functions  can  then  be  used  either  for  direct  analysis  or  to 
recover  a  predictive  distribution.  We  give  an  overview  of  Beta  processes  and, 
more  generally,  the  stochastic  process  approach  to  Bayesian  nonparame tries. 
We  also  include  a  practical  description  of  how  Beta  processes  can  be  simulated. 

Simulation  in  a  Bayesian  setting  requires  accurate  generators  that  can 
produce  variates  from  standard  and  truncated  probability  distributions  for  use 
in  specialized  algorithms.  We  have  therefore  written  a  C++  class  library  that 
will  generate  random  variates  for  use  in  a  variety  of  contexts.  We  needed  a 
source  that  is  reliable,  easy  to  use,  and  flexible.  There  are  many  commercial 
packages  available  that  can  produce  quality  random  variates,  but  these  are 
generally  not  flexible  enough  for  use  in  new  (Bayesian)  research  projects.  Our 
C++  class  was  primarily  developed  for  Bayesian  statistical  modeling  but  can 
be  used  anywhere  simulated  random  variates  are  needed. 

Included  is  a  complete  description  of  the  C++  class  usage  and  imple¬ 
mentation.  We  give  specific  examples  and  discuss  each  function  in  detail  with 
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general  explanations  and  references  to  the  underlying  algorithm.  The  code  has 
been  thoroughly  tested  and  debugged,  but  we  assume  that  the  user  has  a  basic 
understanding  of  probability  theory.  Given  the  usage  description,  references, 
and  attached  source  code,  one  should  be  able  to  incorporate  these  functions 
into  almost  any  simulation  study. 

The  report  is  organized  as  follows:  in  chapter  2  we  discuss  random 
variate  generation  and  the  implementation  and  usage  of  our  C++  class  library. 
Chapter  3  includes  the  theory  and  computational  algorithms  associated  with 
some  Bayesian  nonparametric  methods,  and  chapter  4  concludes  the  report 
with  our  summary  and  possible  areas  for  future  research.  The  appendices 
include  source  code  from  our  C++  class  library. 
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Chapter  2 


Pseudo-Random  Variate  Generation  in  CH — \- 


When  performing  Bayesian  statistical  analysis,  nonparametric,  para¬ 
metric,  and  mixed,  one  frequently  rnnst  simulate  draws  from  standard  proba¬ 
bility  distributions,  predictive  distributions,  and/or  stochastic  processes.  The 
need  arises  in  building  simulation  models  and  while  using  sampling  techniques 
when  non-conjugate  set-ups  render  posterior  distributions  intractable.  In  ei¬ 
ther  case  one  needs  a  good  source  of  uniform  and  non-uniform  pseudo-random 
variates  from  known  distributions.  We  use  the  term  “pseudo”  here  as  true  ran¬ 
domness  is  impossible  from  a  computer;  generated  data  that  appears  random 
is  sufficient  when  trying  to  glean  information  in  a  stochastic  setting. 

We  wanted  a  random  variate  (RV)  source  that  is  reusable,  adaptable, 
and  reliable  so  that  it  can  be  used  in  future  specialized  Bayesian  simulation 
projects.  Thus  we  built  a  class  library  in  C++  to  handle  all  random  num¬ 
ber  generation,  some  statistical  processing  of  data,  and  specialized  Bayesian 
functions  to  include  variate  generation  from  truncated  densities  using  auxil¬ 
iary  variable  techniques  [see  Damien  et  ah,  1999].  C++  was  chosen  for  its 
versatility  and  object  oriented  nature.  For  example,  all  statistical  functions, 
as  well  as  multivariate  generator  functions,  can  be  passed  vectors  from  the 
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C++  standard  library  “vector”  class.  Using  vectors  in  lien  of  arrays  results 
in  code  that  is  faster/more  efficient  and  does  not  require  the  user  to  pass  vec¬ 
tor  sizes,  which  makes  functions  more  general.  In  this  chapter  we  discuss  the 
implementation  and  usage  of  functions  from  our  random  variate  class  library. 

Section  2.1  contains  a  general  description  of  how  to  implement  the  class 
and  its  member  functions.  This  includes  detailed  sample  code  that  should 
make  the  usage  easier  to  understand.  Sections  2.2  and  2.3  discuss  the  class 
member  functions  associated  with  uniform  and  non-uniform  variate  generation 
respectively.  In  section  2.4  we  discuss  algorithms  and  associated  functions  for 
sampling  from  truncated  densities.  We  conclude  the  chapter  in  section  2.5 
with  the  statistical  data  handling  functions. 

2.1  Class  Implementation  and  Usage 

Our  random  variate  class  library  consists  of  four  hies:  RanV.h,  RanV.cpp, 
mrand_seeds.h,  and  BayesRV.cpp.  In  the  usual  manner  RanV.h  contains  the 
class  declarations,  and  RanV.cpp  and  BayesRV.cpp  contain  the  associated 
source  code.  The  header  hie,  mrand_seeds.h,  contains  the  10,000  six  vector 
seeds  needed  for  the  random  number  generators.  To  implement,  place  the 
header  hies  (*.h)  in  an  appropriate  directory  so  that  the  compiler  can  locate 
them,  and  compile  and  link  RanV.cpp  and  BayesRV.cpp  with  the  .cpp  hie 
containing  function  main(). 

Once  RanV.h  is  included  in  the  main  .cpp  hie  with  the  statement 
#include  ‘  'RanV.h’  ’,  merely  instantiate  objects  of  type  RanV  or  BayesRV 
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(with  the  statements  RanV  obj ect  ()  or  BayesRV  obj ect  () )  to  use  the  classes’ 
member  functions.  To  set  a  particular  seed  value,  instantiate  objects  with  in¬ 
put  in  the  interval  [1,10000].  The  objects’  constructors  will  set  the  appropriate 
seed  value.  At  any  other  point  the  seed  can  be  set  to  a  new  value  with  a  call 
to  object .  setSeed(int)1.  Seed  will  default  to  12  otherwise.  C++  source 
code  from  RanV.cpp  is  attached  in  appendix  A  and  code  from  BayesRV. cpp 
in  appendix  B.  Header  hie  RanV.h  is  in  appendix  C. 

The  following  is  complete  sample  C++  program  that  shows  exactly 
how  classes  RanV  and  BayesRV  are  implemented.  The  program  instantiates 
objects  of  type  RanV  and  BayesRV,  calls  functions  from  each  class,  and  outputs 
results  to  a  computer  monitor. 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 


1We  assume  “object”  to  be  the  general  object  instantiated  with  the  statements  RanV 
object  or  BayesRV  object. 


//Sampl.cpp,  Sample  program  that  implements 
//classes  RanV  and  BayesRV 

#include<iostream> 

#include<vector> 

using  namespace  std; 

#include  "RanV.h" 

int  mainO 
{ 

//declare  variables: 
double  U,  N,  Chi,  Gam,  TRe; 
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16 


RanV  rv(55)  ; 


//instantiate  RanV  object 


17 

is  BayesRV  bay(155) ;  //instantiate  BayesRV  object 

19 

20  //call  class  member  functions 


21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 


U  =  rv.mrandO  ; 

N  =  rv .  stdnormO  ; 

Chi  =  rv . ChiSquare (8) ; 
Gam  =  rv . gamma (2 . 4) ; 


//uniform (0, 1)  variate 
//normal (0 , 1)  variate 

//ChiSquare  n  =  8  d.f. 
//Gamma (a),  a  =  2.4 


//return  truncated  exponential  with  rate  =3.0 
//and  truncation  limits  [1,3] 

TRe  =  bay . TRexpn(3 . 0, 1 , 3) ; 


//output  results 

cout  <<  U  <<  endl  <<  N  <<  endl  <<  Chi  <<  endl  <<  Gam 
<<  endl  <<  endl 
<<  TRe  <<  endl; 


35 

36  return  0;  //successful  termination  of  program 

37  }//end  function  main 


In  lines  4-5  we  include  header  files  from  the  standard  template  library. 
We  need  “iostream”  since  we  are  using  cout  and  endl,  and  “vector”  is  neces¬ 
sary  since  class  RanV  employs  vectors.  Line  7  is  necessary  since  cout  and  endl 
are  both  from  the  “std”  namespace.  Line  9  contains  the  #include  statement 
for  the  RanV  header  file,  and  line  11  is  the  official  start  of  the  program. 

The  program  itself  is  well  commented  (lines  beginning  with  “/ /”).  The 
only  items  of  note  are  in  lines  16  and  18.  These  are  where  our  RanV  and 
BayesRV  objects  are  instantiated.  Note  here  that  the  seeds  for  the  uniform 
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generator  for  each  object  are  set  to  55  and  155  respectively.  At  any  subse¬ 
quent  point  during  the  program,  one  could  set  the  seeds  to  a  new  value  with 
rv . setSeed(int)  or  bay . setSeed2(int) . 

2.2  Uniform  Variate  Functions 

All  non-uniform  random  variate  (RV)  generators  in  this  class  library 
are  implemented  using  various  methods  (e.g.  transformations,  acceptance- 
rejection)  involving  Uniform(0,l)  RVs.  It  was  therefore  imperative  that  we 
start  by  developing  a  good  source  of  uniform  variates.  “Good”  uniform  variates 
are  those  that  are  evenly  distributed  across  (0,1)  in  all  dimensions  and  appear 
to  be  independent. 

For  Uniform(0,l)  random  numbers  we  employed  a  composite  Multiple 
Recursive  Generator  (MRG)  as  presented  in  Law  and  Kelton  [2000]  and  de¬ 
veloped  by  L’Ecuyer  [1999].  The  uniform  RV  source  combines  2  MRGs  and  is 
defined  by: 

ZM  =  (1, 403,  580  ZM_ 2  -  810,  728  ZM_3)[mod(232  -  209)] 

Z2,i  =  (527,  612  Z2li-i  -  1,  370,  589  Z2>i_3)[mod( 232  -  22,  853)] 

Vi  =  {Zi,i  -  Z2ti)[mod{ 232  -  209)] 

U=  Y i 

1  232  -  209 

The  parameters  were  chosen  carefully  by  L’Ecuyer  [1999],  and  the  generator 
has  period  of  approximately  2191  with,  according  to  Law  and  Kelton  [2000], 
good  statistical  properties  through  dimension  32.  The  resulting  variates  passed 


all  statistical  tests  for  uniformity  and  independence  for  which  we  subjected 
them.  The  source  code  includes  a  header  hie  with  10,000  6-vector  seeds  spaced 
1016  apart  which  makes  the  generator  a  prime  candidate  for  parallel  computing 
(although  parallel  processors  were  not  used  in  this  application). 

There  are  four  functions  associated  with  Uniform  variate  generation: 
mrandO  ,  unif(),  mrandstO,  and  mrandgt  ().  When  calling  object  .mrandO 
no  input  is  required;  the  seed  value  is  set  when  '  'object’  ’  is  instantiated,  or 
with  the  function  setSeed (stream)  where  “stream”  is  an  integer  on  [1,10000]. 
The  output  will  be  a  uniform  RV  on  the  interval  (0,1).  Figure  2.1  shows  a 
histogram  of  64,000  variates  from  mrandO2.  Simply  use  object  .unit  (a,  b) 
with  real  valued  (a,b)  to  obtain  general  uniform  RVs  on  the  interval  (a,b). 


Figure  2.1:  Uniform(0,l)  random  variates 

Functions  mrandstO  and  mrandgt ()  can  be  used  to  set  and  get  re- 

2The  horizontal  axis  contains  partitioned  X  values,  and  the  vertical  axis  contains  the 
respective  frequencies  of  the  sample. 
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spectively  the  actual  6-vector  seed  used  in  the  uniform  generator.  Call  the 
function  object  .mrandst  (vector<double>  seed,  str)  with  a  vector  con¬ 
taining  6  real  numbers  and  an  integer  on  [0,10000]  to  set  6-vector  stream  “str” 
to  that  contained  in  vector  “seed” .  To  get  the  latest  6- vector  used  in  the  gen¬ 
erator,  pass  an  empty  vector  and  the  stream  number  using  object  ,mr an dgt 
(vector<double>  seed,  str). 

2.3  Non-Uniform  Variate  Functions 

Given  a  good  source  of  Uniform(0,l)  variates,  most  common  non-uniform 
RVs  are  easily  generated.  Published  research  in  this  area  is  plentiful,  robust, 
and  well  developed.  For  an  excellent,  comprehensive  source  one  can  turn  to 
Devroye’s  Non-uniform  Random  Variate  Generation  [see  Devroye,  1986].  The 
book  is  now  out  of  print,  but  an  Adobe  pdf  version  has  been  posted  by  the  au¬ 
thor  and  can  be  found  at  http://jeff.cs.mcgill.ca/~luc/rnbookindex. 
html. 

Our  class  library  can  be  implemented  to  generate  data  from  the  follow¬ 
ing  well-known  continuous  and  discrete  distributions:  exponential,  Weibull, 
gamma,  beta,  Dirichlet,  normal,  lognormal,  Chi-square,  Bernoulli,  binomial, 
multinomial,  and  Poisson.  The  only  continuous  distribution  that  is  imple¬ 
mented  via  direct  inversion  is  the  exponential.  This  comes  from  the  result 
that  Exponential (cc)  variates  can  be  generated  by  X  =  —  ^  ln(U )  where  U  is 

Uniform(0,l).  Sums  of  exponentials  can  then  be  used  to  generate  Gamma(a) 

1 

where  a  €  N  and  Weibull  RVs  are  simply  X&  where  X  is  an  exponential  RV. 
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The  respective  nonuniform  probability  density  and  mass  functions  (pdf’s 
and  prnf’s)  implemented  in  this  package  are: 

1.  Exponential  A):  fxix)  —  \e~Xx,  x  >  0 

2.  Weibull(a,  A):  /,v(x)  =  aAai“_1  x  >  0 

3.  Gamma(a):  fx(x)  =  yyyy  xa~2  e~x,  x  >  0 

4.  Beta(a,  (3):  fx(x)  =  xa~l  (1  -  x)^1,  0  <  x  <  1 

5.  Dirichlct(ai, . . .  ,ak):  fx1,...,xk(x1, . , .  ,xk)  o J2xi  =  l 

-j  (a;  — /z)^ 

6.  Normal (/i,  cr2):  fx(x)  =  e  2ct2  ,  —  oo  <  x  <  oo 

-i  (ln(a:)  — 

7.  Lognormal  (/i,  cr2):  fx(x)  =  ,  x>0 

8.  Chi-square (r):  fx(x)  =  r(r/2)12(r/2)  ar^/2)-1  e~a/2,  a;  >  0 

9.  Bernoulli  (p):  fx{x)  =  px  (1  —  p)1_a:,  x  =  0, 1 

10.  Binomial (n, p):  /x(z)  =  yr^zyyr  px  (1  -  p)1-*,  x  =  0, 1, . . . , n 

11.  Multinomial  (pi, . .  .  ,pfc):  fx1,...,xk(x1,  ...,xk)<x  pT  ■  ■  ■  pxk  ,  &*  =  1 

12.  Poisson(A):  /v(a;)  =  e"A^,  x  =  0, 1, . . . 

To  obtain  exponential  variates,  use  the  function  expn (lambda) 3  with 
any  lambda  >  0.  Weibull  variates  are  obtained  with  the  function  weibull  (shape , 

3From  this  point  forward,  we  will  assume  the  reader  knows  to  call  a  function  using 
“object”  and  drop  the  object .  functionO  notation. 
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scale)  where  the  shape  and  scale  parameters  are  both  real  valued  and  greater 
than  0.  Figure  2.2  contains  a  histogram  of  64,000  exponential(2.0)  variates  and 
figure  2.3  a  plot  of  Weibull(2. 0,4.0)  variates4. 


Figure  2.2:  Exponential  variates:  rate  =  2.0 

Standard  normal  variates  (Normal(0,l))  are  generated  using  the  polar 
method  [see  Law  and  Kelton,  2000]  and  are  easily  converted  to  the  general 
case,  Normal (p,  cr2),  via  Y  —  X  *  a  +  /i.  Standard  normals  can  also  be  used  to 
generate  lognormal  and  chi-square  variates  through  Y  =  ex  and  Y  =  J^X2 
respectively. 

Function  stdnormO  requires  no  input  and  returns  a  Normal(0,l)  ran¬ 
dom  variate.  See  figure  2.4  for  a  histogram  of  generated  standard  normal  vari¬ 
ates.  Normal(p,  a2)  random  variates  are  generated  using  norm  (mu,  sig2), 

4A11  histogram  plots  presented  in  this  chapter  use  64,000  generated  values.  Also,  hori¬ 
zontal  axes  contain  partitioned  X  values,  and  vertical  axes  contain  respective  frequencies  of 
the  sample. 


12 


Figure  2.3:  Weibull  variates:  shape  =  2.0,  scale  =  4.0 

and  lognormal  variates  are  generated  using  function  logn(mu,  sig2).  Pa¬ 
rameters  mu  and  sig2  must  be  real  valued  with  sig2  >  0  in  both  cases  and 
mu  >  0  for  lognormal.  To  obtain  a  Chi  Squared  random  variate  with  n  degrees 
of  freedom,  call  function  ChiSquare(n)  with  integer  valued  n  >  1. 

For  the  generalized  Gamma(a),  a  e  R,  a  >  0,  we  use  a  fast  acceptance- 
rejection  routine5.  For  a  >  1  we  use  a  routine  developed  by  Cheng  [1977]  and 
presented  in  Law  and  Kelton  [2000].  For  a  <  1  we  use  a  routine  attributed 
to  Ahrens  and  Dieter  (1974)  [see  Law  and  Kelton,  2000].  Gamma(a)  RVs  can 
then  be  used  to  obtain  Beta(au,  0:2)  through  Y  =  X^X2  and  Dirichlet  through 
Y  =  Gamma(a,  f3)  is  simply  (5  X ,  where  X  =  Gamma(a).  See  figure  2.5 
for  a  histogram  plot  of  Gamma(5.4)  variates. 

5Note  that  in  the  case  of  very  small  a,  no  acceptance-rejection  routine  will  be  fast.  We 
cover  an  algorithm  for  this  special  case  in  section  2.4. 
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Figure  2.4:  Standard  Normal  Variates:  mean  =  0,  variance  =  1 

A  Gamma(a)  variate  is  returned  with  the  function  gamma  (alpha)  where 
input  alpha  must  be  real  valued  and  greater  than  0.  Again,  to  obtain  a 
Gamma(a,  j3)  variate,  multiply  a  Gamma(a)  by  (5.  Beta  variates  are  ob¬ 
tained  using  function  beta(al,  a2)  with  real  valued  al  >  0  and  a2  >  0.6 
Now,  the  Dirichlet  distribution  is  a  multivariate  generalization  of  the  Beta 
distribution.  For  a  vector  of  realizations  from  a  Dirichlet  distribution  you 
must  pass  a  parameter  vector  and  an  empty  vector  for  x  values  with  the 
function  Dirichlet  (vector<double>  P,  vector<double>  X).  The  parame¬ 
ters  in  vector  P  must  sum  to  1.  After  the  function  is  called,  vector  X  will 
contain  the  desired  random  vector. 

The  discrete  distributions  that  we  implemented  here  are  Bernoulli,  bi¬ 
nomial,  multinomial,  and  Poisson.  Bernoulli  RVs  are  simple  transformations 

6If  al  and/or  a2  are  very  small,  use  the  beta  function  described  in  section  2.4. 
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Figure  2.5:  Standard  Gamma  Variates:  alpha  =  5.4,  beta  =  1 

of  uniforms,  binomials  are  sums  of  Bernoullis,  and  multinomial  random  vec¬ 
tors  are  obtained  using  binomials.  Poisson  RVs  are  a  little  bit  tricky.  We 
employ  an  algorithm  developed  by  Press  et  ah  [2002]  that  separates  cases  into 
A  <  12.0  and  A  >=  12.0.  For  A  <  12.0  the  code  directly  exploits  the  Poisson 
distribution’s  relationship  with  the  exponential.  When  A  >=  12.0  we  use  an 
acceptance-rejection  routine  to  control  processing  time  as  A  gets  large.  See 
figure  2.6  for  a  plot  of  generated  Poisson  variates. 

Function  bern(p)  is  used  to  generate  a  Bernoulli  variate  with  param¬ 
eter  0  <  p  <  1.  To  generate  a  Binomial  variate  call  function  binom(n, 
p)  with  integer  input  n  >  2  and,  again,  0  <  p  <  1.  Now,  the  multino¬ 
mial  function  is  multivariate  and  also  requires  two  vectors  for  input:  one  for 
parameters  such  that  all  sum  to  1  and  one  for  return  of  x  values.  Calling 
multnom(vector<double>  P,  vector<int>  X)  results  in  the  desired  random 
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Figure  2.6:  Poisson  Random  Variates:  mean  =  8.4 

vector  being  stored  in  vector  X.  Finally,  to  obtain  Poisson  variates  (as  in  fig¬ 
ure  2.6),  use  the  function  Poisson(lambda)  with  lambda  >  0. 


2.4  Sampling  from  Truncated  Densities 

In  Bayesian  simulation  one  frequently  must  sample  from  truncated  den¬ 
sities  and/or  distributions  with  exceptionally  narrow  spread.  In  both  cases, 
using  an  acceptance-rejection  routine  is  not  practical.  In  fact,  if  the  trunca¬ 
tion  interval  and/or  spread  is  very  small,  obtaining  a  sample  using  acceptance- 
rejection  may  be  impossible.  To  sample  from  such  distributions  we  have  im¬ 
plemented  in  C++  algorithms  employing  the  auxiliary  variable  technique  de¬ 
veloped  by  Damien  et  al.  [1999].  In  this  section  we  describe  three  algorithms 
from  Damien  and  Walker  [2001]  for  sampling  truncated  Beta,  Gamma,  and 
Bivariate  Normal  densities. 
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The  auxiliary  variable  technique  was  developed  for  ease  of  use  in  a 
Gibbs  sampler  [see  Casella  and  George,  1992],  In  general,  latent  (auxiliary) 
variables  can  be  strategically  added  to  probability  distributions  so  that  the 
full  conditionals  are  known  and  easy  to  sample.  The  latent  variables  must 
be  added  such  that  the  marginal  densities  include  the  desired  distribution. 
In  this  manner  a  Gibbs  sampler  is  used  with  the  full  conditionals  to  produce 
variate  (s)  from  the  original  distribution. 

Introducing  one  auxiliary  variable  to  the  Beta  distribution  reduces  the 
complicated  task  of  generating  truncated  Beta  variates  to  merely  sampling 
a  Uniform  and  an  exponential  function.  The  truncated  Beta  distribution  is 
given,  up  to  proportionality,  by 

fx{x)  oc  xa~l  (1  —  x)13 /(i6  (a,  b )) 

where  a,  /3  >  0  and  0  <  a  <  b  <  1  and  /  is  the  indicator  function.  Damien 
and  Walker  [2001]  introduce  latent  variable  Y  such  that  the  joint  density  is 
given  by 

fx,r(x,y )  oc  xa~l  I  (y  <  (1  -  xf^1,  x  G  (a,  b)).  (2.1) 

Note  that 

(1  —  x)P— 1 

a;0”1 1(y  <  (1  -  a;)^1)  dy  =  a;""1  (1  -  xf~\ 

i.e.  the  marginal  density  of  X  is  again  the  Beta  distribution. 

Given  the  joint  density  in  (2.1),  the  full  conditional  for  Y  given  X  is 
simply  Uniform(0,  (1  —  a:)'3-1).  The  full  conditional  for  X  given  Y  depends 
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on  whether  /?  <  1  or  /?  >  1;  if  /?  =  1  the  Beta  distribution  can  be  sampled 
directly  using  inversion.  So  for  (3  >  1  the  full  conditional  of  X  is 

fx\v(x\y)  oc  xa^l{x  G  (a,min{6, 1  —  2/1/(/3_1)}) )  -  (2.2) 

Now  when  (3  <  1  the  full  conditional  of  X  becomes 

fx\y(x\y)  oc  xa~l  l(x  G  (max{a,  1  -  ?/1/(/3_1)},  b )).  (2.3) 

Both  can  be  sampled  directly  using  the  inverse  CDF  technique  with 

F{x)  =  - - —I(xe{a,b)). 

w  ba  —  aa  v  v  ” 

Now,  the  Gibbs  sampler  algorithm  for  generating  a  truncated  Beta 
reduces  to: 

1.  INITIALIZE  X 

2.  FOR  5  ITERATIONS 

GENERATE  Y  AS  UNIFORM(0,  (1  -  a;)^-1) 

GENERATE  X  FROM  (2.2)  OR  (2.3) 

3.  RETURN  X 

The  biggest  advantage  of  this  algorithm  over  acceptance-rejection  tech¬ 
niques,  if  there  are  no  truncation  limits,  occurs  when  a  and/or  (3  are  very  small. 
See  figure  2.7  for  a  sample  of  Beta(0.2,0.2)  variates.  These  are  generated  using 
function  TRbeta (alpha,  beta,  a,  b)  from  the  “BayesRV”  class  where  a  and 
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Figure  2.7:  Beta  Random  Variates:  alpha  =  0.2,  beta  =  0.2 

b  are  the  truncation  limits.  To  obtain  non  truncated  variates,  merely  set  a  = 
0  and  b  =  1. 

The  method  for  truncated  Gamma  variates  is  very  similar  to  that  of 
Beta.  The  Gibbs  sampler  algorithm  is  exactly  the  same  as  above  but  with 
slightly  different  conditionals.  The  truncated  standard  (/3  =  1)  Gamma  den¬ 
sity  is  given  to  proportionality  by: 

fx(x)  oc  x a_1  exp(— x)  l(x  G  (a,  &)) 

where  0  <  a  <  b  <  oo.  Damien  and  Walker  [2001]  introduce  latent  variable  Y 
such  that  the  joint  density  is  given  as: 

Ix,y(x,  y)  OC  xa_1 1 {%/  <  exp(— x),  x  G  (a,  6)).  (2.4) 

Again,  note  that,  integrating  with  respect  to  y: 

rexp(—x) 

/  x“_1  l(y  <  exp(-x))  dy  =  x“_1  exp(— x). 

Jo 
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Now,  given  the  joint  density  in  (2.4)  we  have  the  full  conditional  for  Y 
given  A"  as  Uniform(0,  exp(— x)).  The  full  conditional  for  X  given  Y  is  given 
as 

fx\r(x\y)  oc  xa~l  l{x  e  (a,mm{b,-log(y)})).  (2.5) 

Equation  (2.5)  is  very  similar  to  that  of  (2.2)  and  (2.3)  and  can  be  sampled 
using  inversion. 

To  obtain  a  truncated  Gamma  variate  call  function  TRgammaC alpha, 
a,  b)  from  the  “BayesRV”  class.  To  obtain  non  truncated  variates  simply 
call  TRgammaC )  with  a  =  0  and  b  arbitrarily  high  compared  to  alpha.  Again, 
the  advantages  over  the  standard  Gamma  function  occur  when  alpha  is  very 
small.  See  figure  2.8  for  a  sample  of  Gamma(0.57)  variates. 

Like  the  truncated  Beta  and  Gamma  distributions,  sampling  from  a 
truncated  multivariate  Normal  density  can  be  greatly  simplified  using  the  aux¬ 
iliary  variable  technique  of  Damien  et  al.  [1999].  Without  this  technique  the 
task  can  be  arduous  at  best.  Here  we  have  employed  an  algorithm  from  Damien 
and  Walker  [2001]  to  handle  the  bivariate  case.  Given  a  good  matrix  handling 
C++  library,  this  method  can  easily  be  extended  to  larger  random  vectors. 

The  truncated  multivariate  Normal  density  is  given  to  proportionality 
by 

fxu...,xp(xi, . .  .,xp)  oc  exp(  -  -  (x-  +)'E“1  (x  -  /z))  l(x  e  A) 
where  /i  is  the  mean  vector  and  E  is  the  covariance  matrix.  Damien  and  Walker 
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[2001]  introduce  only  one  latent  variable  Y  and  define  the  joint  density  as 

y 

fxu...,xp,Y{x1,..  .,xp,y)  oc  exp(  -  -)  I  (y  >  (x  -  /i)'  E_1  (x-/i),i£  A) 
where  the  full  conditional  for  each  Xt  is  given  as 

fxi\x-i,Y(.Xi\x-i,y)  oc  I(xiEAi). 

The  conditionals  are  then  uniform  densities  with  Ai  =  (a*,  A)  D  B,  where  Bi  is 
the  set  {xi\x-i  :  (x  —  (i)'  £-1  (x  —  ft)  <  y}.  The  bounds  for  Bi  are  therefore 
found  by  solving  a  quadratic  equation.  Also,  the  full  conditional  for  Y  given 
X  is  a  truncated  exponential  distribution. 


Figure  2.8:  Gamma  Random  Variates:  alpha  =  0.57,  beta  =  1 

Since  the  truncated  exponential  distribution  is  used  frequently  in  sim¬ 
ulation  algorithms,  other  than  just  sampling  truncated  multivariate  normal, 
we  created  the  function  TRexpn( lambda,  a,  b).  Call  this  function  with  real 
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b  =  3.0 

valued  lambda  greater  than  0  and  0<a<b<ooto  obtain  truncated  expo¬ 
nential  variates  as  in  figure  2.9. 

Now,  to  obtain  truncated  bivariate  Normal  variates  call  the  function 
TRbiNorm(X,  Mu,  Sigma,  A).  X  is  an  empty  vector<double>  to  hold  the 
returned  variates,  Mu  is  a  vector<double>  containing  the  means,  and  Sigma 
is  a  2x2  double  subscripted  array  containing  the  covariance  matrix.  A  is  again 
a  2x2  double  subscripted  array.  The  first  row  contains  the  truncation  limits 
for  Xi  and  the  second  row  contains  the  truncation  limits  for  A^. 

Using  function  TRbiNormO  we  reproduce  a  simulation  described  in 
Damien  and  Walker  [2001]  of  10,000  samples  from  a  bivariate  Normal  den¬ 
sity  truncated  on  the  unit  circle  centered  at  (0.5, 0.5).  The  means  here  are  0 
and  E  is  given  with  unit  variances  and  covariance  =  0.9.  Figure  2.10  contains 
a  graph  of  the  produced  variates.  Notice  the  truncation  limits  and  the  positive 
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correlation  between  X\  and  X2. 


Figure  2.10:  10,000  Truncated  Bivariate  Normal  Variates 

2.5  Statistical  Data  Handling  in  C-| — b 

In  addition  to  pseudo-random  variate  generation,  we  can  use  our  class 
library  to  input  vectors  of  data  from  sequential,  delimited  data  files  and  per¬ 
form  empirical  statistical  analysis.  Analysis  includes  calculating  population 
mean  and  variance,  covariance  between  two  vectors  of  data,  and  empirical 
quartiles.  Input  is  stored  in  vectors  so  all  data  handling  functions  can  be  used 
with  internally  simulated  values  as  well. 

External  data  hies  can  be  delimited  in  almost  any  manner  to  include 
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spaces,  commas,  lines,  and  tabs.  To  get  external  data,  call  the  function 
getdata(vector<>  D,  filename)  with  an  empty  vector  and  filename  to  in¬ 
clude  the  hie  extension,  i.e.  “hlename.txt”  or  “hlename.dat”.  The  data  hie 
must  be  in  the  same  directory  as  the  compiled  executable.  The  data  will  be 
stored  in  vector  “D”. 

Data  handling  functions  include  cmom  (vectoro  D ,  mean ,  var ) ,  cov  ( 
vectoro  Dl,  vectoro  D2,  cov) ,  and  quant  (vectoro  D,  vectoro  Q). 
The  inputs  to  cmom()  are  a  data  vector,  D,  and  double  valued  variables  mean 
and  var  to  hold  the  outputs  for  population  mean  and  population  variance 
respectively.  Function  cov()  requires  two  data  vectors,  Dl  and  D2,  and  double 
valued  cov  to  hold  the  output  sample  covariance  between  the  two  data  vectors. 

Finally,  function  quantO  requires  a  data  vector,  D,  and  a  vector  of 
quantiles,  Q.  For  output,  the  data  vector  is  returned  numerically  sorted,  and 
function  quant  ()  replaces  the  percentile  values  in  the  quantile  vector  with  the 
desired  data  value.  For  example,  Q  might  contain  as  input  {0 . 1 , 0 . 5 , 0 . 9} 
to  return  the  10th,  50th,  and  90th  percentiles  respectively.  If  the  data  vec¬ 
tor  contained  standard  Normal  variates,  the  quantile  vector  above  would  be 
returned  as  approximately  {-1 . 282 ,0,1. 282}. 
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Chapter  3 


Bayesian  Nonparametrics 

3.1  Introduction 

The  focus  of  the  Bayes  portion  of  this  report  is  on  the  use  of  Bayesian 
nonparametric  techniques  for  reliability  (lifetime  data)  analysis.  In  general  we 
would  like  to  glean  useful  information  and  make  accurate  predictions  using 
lifetime  data  collected  about  an  observed  set  of  items.  In  a  traditional  para¬ 
metric  setting  we  would  make  certain  assumptions  about  the  family  of  dis¬ 
tributions  from  which  these  lifetimes  arise  and  use  collected  data  to  estimate 
the  unknown,  but  fixed,  parameter(s)  from  the  chosen  distribution  family.  In 
Bayesian  parametric  analysis  we  view  the  parameter(s)  as  random  and  assign 
prior  distribution (s)  from  which  we  make  inferences  about  the  parameter  (s) 
in  question.  We  then  use  Bayes  theorem  to  make  data  conditional  updates  to 
the  prior  distribution(s). 

Traditional  nonparametric  techniques  [see  Siegel  and  Castellan  Jr.,  1988] 
involve  not  making  any  assumptions  about  the  family  of  distributions  and  us¬ 
ing  the  data  to  build  unknown,  yet  fixed,  distributions  from  which  to  make 
predictions.  Now,  in  parallel  with  the  parametric  setting,  Bayesian  nonpara¬ 
metric  methods  “randomize”  the  distribution  in  question  and  place  a  prior  to 
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which  we  can  make  data  conditional  updates  (posterior).  We  can  then  use  the 
posterior  to  make  inferences  about  the  distribution  and  in  turn  make  predic¬ 
tions  about  lifetimes  in  question.  A  current  survey  of  many  popular  Bayesian 
nonparametric  techniques  with  follow-on  discussions  can  be  found  in  Walker 
et  al.  [1999]. 

We’ve  chosen  to  investigate  the  use  of  stochastic  process  priors  in 
Bayesian  nonparametric  models.  Since  Ferguson  [see  Ferguson,  1973]  intro¬ 
duced  the  Dirichlet  process  prior  in  his  seminal  1973  paper,  the  stochastic 
process  approach  to  Bayesian  nonparametrics  has  become  very  popular  for  its 
stable  theoretical  grounding  and  ease  of  use  [see  Walker  et  ah,  1999].  Ferguson 
[1973]  showed  that  the  Dirichlet  process  has  large  support  and  is  conjugate  to 
the  space  of  cumulative  distribution  functions  (CDFs),  i.e.  the  data  condi¬ 
tional  posterior  process  is  again  a  Dirichlet  process. 

In  the  stochastic  process  approach  to  Bayesian  nonparametrics,  we  de¬ 
fine  a  stochastic  process  with  index  set  on  the  sample  space,  fl,  where  spe¬ 
cific  realizations  are  cumulative  distribution  functions.  In  this  setting,  if  the 
stochastic  process  is  carefully  chosen  (as  in  the  Dirichlet  process  introduced 
above),  updates  can  be  made  using  observations  from  the  distribution  in  ques¬ 
tion  that  will  result  in  a  well  defined  posterior  process.  The  increments  of  the 
posterior  process  can  then  be  simulated  to  obtain  realizations  of  the  underlying 
CDF,  which  can  in  turn  be  used  to  simulate  predictive  data  sets. 

The  Dirichlet  process  prior  has  one  drawback  in  that  it  leads  to  the 
class  of  discrete  distributions  almost  surely.  To  overcome  this  (while  handling 
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censored  lifetime  data),  Doksum  [1974],  Ferguson  [1974],  and  Ferguson  and 
Phadia  [1979]  developed  a  general  class  of  priors  called  Neutral  to  the  Right 
(NTR)  processes  of  which  the  Dirichlet  is  a  specific  case.  They  showed  that  if 
the  random  distribution  function  F(x),x  E  is  NTR,  then  Z(x)  =  —log(  1  — 
F(x))  is  a  Levy  process  where 

1.  Z(x)  has  non- negative  independent  increments, 

2.  Z(x)  is  non-decreasing  almost  surely, 

3.  Z(x)  is  right  continuous  almost  surely, 

4.  limx^-00[Z(x)\  =  0  almost  surely  and 

5.  limx^+00[Z(x)\  =  oo  almost  surely. 

Ferguson  and  Phadia  [1979]  also  proved  that  if  F(x)  is  NTR,  then  F(x) 
given  a  sample  X i ,  X2 ,  •  •  • ,  Xn ,  possibly  right  censored,  from  F(x)  is  again 
NTR  and  Z(x)  is  again  an  independent  increments  Levy  process.  As  a  result, 
possibly  right  censored  data  can  be  used  to  make  updates  to  the  selected  prior 
process,  then  focus  can  be  turned  to  the  Levy  process  Z ( x )  for  realizations  of 
the  random  distribution  F(x). 

In  this  chapter  we  turn  our  attention  to  a  specific  class  of  Levy  processes 
developed  by  Hjort  [1990],  beta  processes,  which  are  used  as  priors  for  the  space 
of  cumulative  hazard  functions  (CHFs).  The  distribution  function  F(t),t  >  0 
can  be  recovered  from  the  CHF,  A(t),  as  A{t)  =  —log(  1  —  F(t))  in  the  case  of 
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continuous  F(t).  Notice  the  parallel  of  A(t)  to  Z(x)  above.  In  section  3.2  we 
discuss  the  use  of  the  cumulative  hazard  function  in  a  Bayesian  nonparametric 
setting.  Section  3.3  introduces  the  theoretical  framework  behind  the  beta 
process,  and  section  3.4  addresses  the  computational  issues  associated  with 
simulating  the  increments  of  the  beta  process. 

3.2  Random  Cumulative  Hazard  Function 

Most  of  the  emphasis  in  the  Bayesian  nonparametric  literature  has  been 
placed  on  developing  priors  for  the  random  distribution  function  (CDF).  To 
that  effect  the  1970s  saw  the  development  of  Neutral  to  the  Right  (NTR) 
processes  as  priors  for  the  CDF,  F(t),  where  Z(t)  =  —log(l  —  F(t))  is  an  inde¬ 
pendent  increments  Levy  process.  Hjort  [1990],  noticing  a  natural  relationship 
of  hazard  functions  to  Z(t),  chose  to  focus  on  developing  a  stochastic  process 
prior  not  for  the  CDF,  F(t),  but  for  the  cumulative  hazard  function  (CHF), 
A(t),  instead.  According  to  Hjort  [1990]  the  CHF  is  as  basic  to  understand¬ 
ing  survival  phenomenon  as  the  CDF,  and  the  hazard  rate  concept  is  easily 
generalized  to  more  complicated  models. 

We’ve  chosen  only  to  address  the  time  continuous  CHF  model.  Let  T  be 
a  random  variable  with  CDF  F(t)  =  Pr(T  <  t )  on  [0,  oo)  and  F(0)  =  0.  The 
cumulative  hazard  rate  for  T  is  a  nonnegative,  nondecreasing,  right  continuous 
function  A  on  [0,  oo)  where 

dA(s)  =  A[s, s  +  ds)  =  Pr{T  e  [s, s  +  ds)\T  >  s}  = 
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dF(s) 
F[s,  oo) 


Now  Hjort  defines  for  0  <  a  <  b  <  oo 


A[a,  b ) 


M) 


dF(s) 

F[s,  oo) 


and  requires 

F[a,b)  —  /  F[s,  oo)  dA(s). 

J  [a,b) 

We  have  that  A(t)  is  the  value  of  the  function  A  at  the  point  t  such  that  F(t) 
is  recoverable  from  A(t)  by  using  product  integrals[Gil\  and  Johansen,  1990]: 


m  =  i  -  n  {1  —  gL4(s)},  t  >  0.  (3.1) 

«e[o,t] 

If  F  is  continuous  this  reduces  to  the  familiar 


A(t )  =  -log(l  -  F(t)). 


We  now  have  the  framework  for  developing  Hjort’s  prior  distributions 
for  A.  Let  F  be  the  set  of  all  CDFs  F  on  [0,  oo)  having  F(0)  =  0  and  let  B 
be  the  set  of  all  nondecreasing,  right  continuous  functions  B  on  [0,  oo)  having 
B( 0)  =  0.  Define  the  space  of  cumulative  hazard  rates  as 

A  =  {A  G  B  for  which  (3.1)  leads  to  an  F  e  F}. 

To  place  probability  distribution  on  A,  it  is  natural  to  consider  the  non¬ 
negative,  nondecreasing  Levy  processes,  Z(t),  defined  above  where  Z(t)  = 
—log(  1  —  F(t)).  However,  not  every  Levy  process  can  be  used  as  a  prior  for 
A.  To  address  this  limitation  Hjort  constructs  the  beta  process,  whose  sample 
paths  lie  in  A  almost  surely. 
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3.3  The  Beta  Process 


In  this  section  we  outline  the  key  definitions  and  theorems  associated 
with  the  beta  process.  These  include  complicated  Levy  representations  that 
define  the  increments  of  the  process.  The  Levy  representation  is  a  formula 
denoting  the  moment  generating  function  (MGF)  or  characteristic  function 
(CF)  of  the  process.  We  will  show  in  the  next  section  how  these  Levy  formulas 
can  be  used  in  practice.  Many  stochastic  processes,  including  most  of  the 
general  processes  developed  as  Bayes  priors  for  random  CDFs  and  CHFs,  have 
no  closed-form  representation  defining  the  distribution  of  the  increments-only 
an  MGF  or  CF.  Hjort’s  beta  process  is  no  exception.  We  should  point  out  that 
this  is  merely  an  overview  of  the  beta  process  theory.  Further  development 
and  proofs  for  all  theorems  can  be  found  in  Hjort  [1990]. 

We  now  define  the  beta  process.  Note  that  the  goal  of  this  definition 
is  to  split  the  process  into  a  continuous  and  a  discrete  component.  This  lends 
itself  to  computational  practicality  as  discussed  in  the  next  section. 

Definition  3.3.1.  Let  Aq  be  a  CHF  with  a  finite  number  of  jumps  taking  place 
at  ti,t2,  •  •  • ,  and  let  c(-)  be  a  piecewise  continuous,  non-negative  function  on 
[0,  cx)).  A  beta  process  with  parameters  c(-),  A(-)  is  a  Levy  process  denoted  by 

A  ~  beta{c(-),  A0(-)},  (3.2) 

if  A  has  Levy  representation 

E[e-eAV]  =  |  ]J  E[e-es-]}  exp{-  (1  -  e-e‘)dMs)},  (3.3) 
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with 


si  =  A{fj}  ~  beta{c(tj )  A0{tj},  c(tj)  (1  -  A0{tj})},  (3.4) 

and 

dLt{s)  =  f  c(z)  s-1  (1  —  dA0:C(z)ds  (3.5) 

Jo 

for  t  >  0  and  0  <  s  <  1,  in  which  A0tC(t)  =  A0(t )  —  ^ ~2t.<t  Ao{tj}  is  A0(t)  with 
its  jumps  removed. 

We  can  restate  the  last  part  of  the  definition  as 

A(t)  =  '^2  Sj  +  Ac(t), 

tj<t 

where  the  jumps,  Sj,  are  independent  beta  random  variables  from  (3.4),  and 
Ac(t)  is  beta{c(-),  A0jC(-)}  as  defined  in  (3.2). 

Assume  that  A0  is  the  prior  guess  at  the  cumulative  hazard,  and  c(s) 
can  be  interpreted  as  the  number  at  risk  at  time  s  in  an  imagined  prior  sample 
with  hazard  rate  corresponding  to  Aq.  The  prior  guess  at  the  cumulative 
hazard  can  be  chosen  from  a  parametric  model  such  as  the  Weibull,  i.e.  Ao(t)  = 
Xa  ta  where  a  is  an  assumed  prior  shape  parameter  and  A  is  an  assumed  prior 
scale  parameter. 

To  define  a  Bayes  posterior  as  developed  by  Hjort  for  the  random  CHF 
A(t),  let  Xi,  X2, . . . ,  Xn  be  a  random  sample  from  F(t)  with  CHF  A  e  A 
(as  defined  above).  Assume  that  (Tj,  <5i), . . . ,  (Tn,  Sn)  are  observed  where  T*  = 
min(Xj,Ci),  <5*  =  I{Xi  <  c,}  and  ci,...,cn  are  censoring  times.  Define  the 
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counting  process  N  and  the  left-continuous  at-risk  process  Y  by 


N{t)  =  I{Ti  <  t  &  <5*  =  1},  Y(t)  =  >  t},  (3.6) 

i= 1  t=  1 

where  I  is  the  indicator  function.  Assume  that  the  censoring  times  are  either 
fixed  or  independent  of  the  lifetimes  Xt.  Also  note  that  dN(t )  =  N(t)  is  the 
number  of  observed  Yds  at  a  particular  time  t.  The  posterior  process  given 
the  data  will  then  be 


Theorem  3.3.1.  (Hjort’s  corollary  j.l)  Let  A  ~  beta{c(-),  A0(-)}  as  in  (3.2). 
Then 


A  |  (Th  Ai), . . . ,  (Tn,  Sn)  ~  betal  c(-)  +  Y(-), 


c{s)  dAo(s)  +  dN(s) 
j  c(s)  +  Y(s) 


(3.7) 


The  posterior  process  contains  fixed  points  of  discontinuity  (at  each  of 
the  observed  times,  censor  and  survival)  even  if  the  prior  does  not. 


3.4  Computational  Issues 

Simulating  sample  paths  from  stochastic  processes  with  no  closed  form 
representation  for  their  increments  can  be  quite  tedious.  The  stochastic  pro¬ 
cesses  described  earlier  in  this  chapter  can  be  separated  into  a  continuous  com¬ 
ponent  and  a  jump  component.  The  jump  components  are  straightforward  in 
most  cases,  as  they  are  usually  drawn  from  standard  probability  distributions. 
The  continuous  component  is  much  more  complicated  in  general.  Damien  et  al. 
[1995]  gives  a  general  method  for  sampling  from  such  distributions  (denoted 
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infinitely  divisible )  and  Damien  et  al.  [1996]  apply  this  method  to  Hjort’s  beta 
process. 

Sampling  from  posterior  beta  processes  can  be  simplified,  though,  by 
approximating  the  continuous  component  with  standard  beta  distributions. 
This  is  a  result  of  the  fact  that  the  accompanying  Levy  measures  Ltfi)  associ¬ 
ated  with  the  beta  process  as  defined  in  (3.2)  are  concentrated  on  the  interval 
(0,1).  The  jump  component  is  already  given  as  a  beta  random  variable  (see 
equation  (3.4)). 

To  simulate  a  beta  process  let  H (f)  denote  a  random  cumulative  hazard 
function  and  partition  the  time  axis  into  times  t  =  0, . .  . ,  T  (a  hirer  grid  results 
in  a  smoother  curve).  We  have  that  H{t )  =  ^ ~2t  <tSj  +  Ac{t)  where  Sj  are 
jumps  at  time  tj  and  Ac[t )  is  the  continuous  component.  Assume  that  the  prior 
process  is  absolutely  continuous,  i.e.  Yhtj<t  —  0  f°r  all  t.  For  the  continuous 
component  we  assume  Ac  ~  Beta(a(t),b(t)),  where  a  =  c{t)  and  b  =  A0(t) 
(the  prior  guess).  Note  that  “Beta”  here  denotes  the  beta  distribution  not  the 
beta  process. 

Take  c(t )  =  k,  where  k:  is  an  integer,  say,  1  or  2  and  A0(t)  =  at.  This 
means  the  underlying  failure  rate  is  linear,  with  mean  A  Choose  a  =  0.1 
for  example.  Now,  to  sample  the  beta  process  PRIOR,  at  each  time  interval 
sample  from  the  appropriate  Beta  (a(t),b(t))  defined  above.  Note  that  the  A0 
parameter  will  change  if  the  time  intervals  are  of  unequal  lengths.  You  can 
also  make  it  more  general  by  setting  c(f)  =  A  t,  where  A  is  a  real  number. 
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Given  data,  the  process  will  have  jumps  at  each  failure  and  censor 
time  tj.  Each  jump  S3  will  have  a  beta  distribution  with  parameters  a  = 
c(t )  A0(tj)  +  dN(t)  and  b  =  c(t )  (1  —  A0(tj)  +  Y(t )  —  dN(t ))  where  dN(t)  is  the 
number  of  failures  by  time  t  and  Y(t)  is  the  number  still  at  risk  by  time  t  as 
defined  in  (3.6). 

Now,  the  a(t)  parameter  in  the  beta  distribution  for  the  continuous 
component  is  updated  by  simply  adding  the  number  of  at  risk  observations 
(right  censored)  from  that  interval  to  c(t).  And  so  a(t)  =  c(t)  +  Y(t).  Up¬ 
dates  to  the  bit)  parameter  of  the  continuous  component  are  slightly  more 
complicated.  First  define  a  couple  of  variables: 

1.  dAo(t)  =  ct;  the  differential  of  A$  given  above.  Note  if  you  change  the 
choice  of  Aq  above,  then  the  differential  will  obviously  change. 

2.  dN(t)  is  the  exact  observations  in  the  interval  of  interest. 

3.  Numerator  =  c(s )  gL40(s)  +  dN(s) 

4.  Denominator  =  c(s)  +  Y (s) 

For  the  first  time  interval  from  time  t  —  0  to  the  first  point  on  the 
grid,  say,  t1,  integrate:  f^0  Nominator  •  This  integrated  quantity  is  your  b(t) 
parameter  in  the  POSTERIOR  beta  process  in  the  first  interval.  To  get  the 
bit)  for  the  second  interval  integrate  from  ti  to  f2,  and  so  on. 
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Key  point:  Note  that  the  a(t)  and  b(t)  parameters  will  likely  change 
from  one  time  interval  to  the  next.  As  you  can  see,  with  complicated  choices 
for  c  and  A0,  the  updates  can  get  nasty  because  the  integral  can  become  awful. 

Using  these  parameters,  simulate  jump  components  and  a  continuous 
component  at  each  time  R  in  the  partition  of  the  time  axis.  The  final  simulated 
value  for  H(t)  will  now  be  ^2tj<t  Sj+Ac(t).  Selecting  the  time  axis  partition  to 
coincide  with  the  EXACT  and  RIGHT  CENSORED  data  points  will  simplify 
things.  To  obtain  a  smoother  estimate  of  the  hazard  function  include  more 
grid  points.  An  estimate  of  the  CDF  F(t)  can  now  be  recovered  using  F{tj)  = 
1  —  [(1  —  F{tj- 1)  *  (1  —  Ac{tj ))  *  (1  —  S(tj))]  as  derived  from  the  product  integrals 
in  (3.1). 
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Chapter  4 


Summary  and  Areas  for  Continued  Research 


Using  Bayesian  nonparametrics  for  simulation  input  modeling  has  ex¬ 
ceptional  potential.  Until  recently  such  methods  were  impractical  due  to  lim¬ 
ited  computing  power  and  relatively  immature  computational  methods.  There 
now  exists  multitudes  of  new  published  research  in  the  area,  and  computing 
power  is  virtually  unlimited.  We  have  given  an  overview  of  Bayesian  non- 
parametric  methods  with  specific  focus  in  hazard  rate  modeling  for  reliability 
analysis.  We  have  also  written  and  presented  two  C++  classes  that  can  be 
used  as  a  foundation  for  almost  any  Bayesian  simulation  project. 

Computer  programs  inherently  are  evolving  entities.  Any  program  or 
algorithm  has  potential  for  improvement.  The  C++  classes  presented  in  chap¬ 
ter  2  are  no  exception.  They  have  been  thoroughly  debugged  and  employ 
reliable,  proven  algorithms  but  are  still  a  work  in  progress.  Also,  they  can 
be  tailored  to  fit  almost  any  project  where  random  variates  are  needed.  The 
classes  are  a  foundation  and  are  meant  to  fill  the  void  between  published 
methodologies  and  actual  working  simulation  programs. 

In  addition,  the  Bayesian  theory  presented  in  chapter  3  presents  a  start¬ 
ing  point  for  future  research.  The  methods  can  be  applied  to  a  “real  world” 
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reliability  problem  and  compared  to  existing  procedures.  Section  3.4  presents 
a  usable  technique  with  examples  for  simulating  Beta  processes,  but  the  algo¬ 
rithm  was  not  actually  coded.  The  examples  should  be  tailored  to  a  specific 
problem  since  prior  knowledge  will  be  different  in  every  case.  The  Bayesian 
nonparametric  method  presented  here  could  be  easily  implemented  in  an  ex¬ 
isting  reliability  study  where  failure  events  are  volatile  and  difficult  to  model. 
Censored  data  is  also  easily  managed  in  this  context. 

For  simulation  input  modeling,  the  hazard  rates  sampled  from  the 
Bayesian  nonparametric  model  can  be  used  to  recover  a  discrete  cumulative 
distribution  function  (CDF)  as  described  in  section  3.4.  This  CDF  can,  in 
turn,  be  used  to  sample  failure  times  for  input  to  a  simulation.  Techniques  for 
sampling  from  a  discrete  CDF  are  well  defined  and  can  be  found  in  Law  and 
Kelton  [2000].  To  obtain  a  smoother  curve  with  more  possible  failure  times, 
simply  create  a  finer  grid  for  sampling  the  hazard  rates. 
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Appendix  A 


Source  Code:  RanV.cpp 


The  following  is  verbatim  source  code  from  RanV.cpp  as  described  in 
Chapter  2. 

//RanV.cpp,  Patrick  J.  Munson,  1  Apr  05 
//Member  function  definitions  for  class  RanV 
//generates  uniform  and  non-uniform  pseudo-random  variates 
//and  performs  some  statistical  data  handling 


#include 

#include 

#include 

#include 

#include 

#include 

#include 


<iostream> 

<cmath> 

<vector> 

<algorithm> 

<f stream> 
"RanV.h" 
"mrand_seeds . h" 


using  namespace  std; 

All  RV  functions  employ  uniform  RVs  from  combined  MRG  as  presented  by 
L’Ecuyer  (1999).  10,000  streams  are  supported,  with  seed  vectors  spaced 

10~16  apart.  Class  RanV  object  must  be  initialized  with  a  seed  value  in 
the  interval  [1,10000],  If  not  seed  defaulted  to  12. 

RanV : : RanV (  int  seed  ) 

{ 

setSeed(  seed  ) ; 

}//end  RanV  constructor 

void  RanV: : setSeed(  int  s  ) 

{ 

stream  =  (  s  >=  1  &&  s  <=  10000  )  ?  s  :  12; 

}//end  function  setSeed 
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//UNIFORM  RVs 

Function  mrandO :  Returns  unit orm(0 , 1) . 

double  RanV : : mrand ( ) 

{ 

//variables 

const  double  ml=4294967087 . 0 ,  m2=4294944443 . 0 ; 
const  double  norm=l . 0/ (ml+1) ; 
int  k; 
double  p; 

double  sl0=drng [stream] [0] ,  sll=drng [stream] [1] ,  sl2=drng [stream] [2] ; 
double  s20=drng [stream] [3] ,  s21=drng [stream] [4] ,  s22=drng [stream] [5] ; 

//calculate  next  "random"  number 

p  =  1403580.0  *  sll  -  810728.0  *  slO; 

k  =  (int) (p/ml) ; 

p  -=  k*ml ;  //p  (mod  ml) 

if  (p  <  0.0)  p  +=  ml; 

slO  =  sll;  sll  =  sl2 ;  sl2  =  p; 

p  =  527612.0  *  s22  -  1370589.0  *  s20; 

k  =  (int) (p/m2) ; 

p  -=  k*m2;  //p  (mod  m2) 

if  (p  <  0.0)  p  +=  m2; 

s20  =  s21;  s21  =  s22;  s22  =  p; 

//update  stream 

drng [stream] [0]  =  slO;  drng [stream] [1]  =  sll;  drng [stream] [2]  =  sl2; 
drng [stream]  [3]  =  s20;  drng [stream] [4]  =  s21;  drng [stream]  [5]  =  s22; 

if  (sl2  <=  s22)  return  ((sl2  -  s22  +  ml)  *  norm); 
else  return  ((sl2  -  s22)  *  norm); 

}//end  function  mrandO 

Function  unifQ:  Returns  uniform(a,b)  . 

double  RanV::unif(  double  a,  double  b  ) 

{ 

double  U; 

U  =  mrandO  ; 

return  U*(b-a)  +  a; 
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}//end  function  unif() 


Function  mrandstO:  Sets  seed  vector  "stream"  to  desired  6-vector. 

void  RanV : : mrandst (  vector<double>  &seed,  int  str  ) 

{ 

for  (  int  i=0;  i<6;  i++  ) 

drng[str][i]  =  seed[i]; 

}//end  function  mrandstO 

Function  mrandgtO:  Returns  the  most  current  6-vector  of  integers. 

void  RanV : : mrandgt (  vector<double>  &seed,  int  str  ) 

{ 

seed.  clearO  ; 

for  (  int  i=0;  i<6;  i++  ) 

seed . push_back(  drng[str][i]  ); 

}//end  function  mrandgtO 
//END  UNIFORM  RVs 


//NONUNIFORM  RVs 

Function  "stdnormO":  Returns  a  N(0,1)  random  variate  using  the 
mrandO  function  and  the  polar  method. 

double  RanV:  :  stdnormO 

{ 

//variables 
static  int  flag=0; 
static  double  X2; 
double  W,  Y,  VI,  V2; 

//avoid  duplicating  effort  using  a  flag,  utilize  both  variates 
if  (  flag  ==  0  )  { 
do  { 

VI  =  2 . 0*mrandO-l .  0 ; 

V2  =  2 . 0*mrandO-l .  0 ; 

W  =  V1*V1  +  V2*V2 ; 

>  while  (W>1.0  II  W  ==  0.0  ) ; 

Y  =  sqrt(  (-2 . 0*log(W) ) /W  ); 

X2  =  V2*Y ; 
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flag  =  1; 
return  V1*Y; 

}  else  -[  //we  have  an  extra  variate  handy — return  it 
flag  =  0; 
return  X2; 

}//end  if/else 
}//end  function  stdnormO 

Function  norm() :  Returns  N(mu,  sig~2) ,  using  stdnormO 

double  RanV::norm(  double  mu,  double  sig2  ) 

{ 

double  nO,  nm; 

nO  =  stdnormO  ; 

nm  =  sqrt (sig2) *n0  +  mu; 

return  nm; 

}//end  function  normO 

Function  lognO :  Returns  lognormal  using  normO  and  a  transformation. 

double  RanV::logn(  double  mu,  double  sig2  ) 

{ 

double  nl,  Nmu,  Nsig2,  lg; 

Nmu  =  log(  mu*mu/ (sqrt (mu*mu+sig2) )  ); 

Nsig2  =  log(  1  +  sig2/(mu*mu)  ); 

nl  =  norm(Nmu,Nsig2) ; 
lg  =  exp(nl) ; 

return  lg; 

}//end  function  lognO 

Function  expnO :  Returns  exponential (lamda)  random  variate  using  mrandO 
and  the  inversion  method. 

double  RanV::expn(  double  lambda  ) 

{ 

double  U,  ed; 
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U  =  mrandO  ; 

ed  =  (-1 . 0/lambda) *log(U) ; 


return  ed; 

}//end  function  expn() 

Function  weibullQ:  Returns  weibull(a,b)  random  variate  using  expn(a) 
and  a  transformation. 

double  RanV: : weibull (  double  shape,  double  scale  ) 
double  a,  w; 
a  =  pow(  scale, shape  ); 
w  =  pow(  expn(a) , (1 . 0/shape)  ); 
return  w; 

}//end  function  weibull () 

Function  bernO  :  Returns  bernoulli  using  mrandO  . 

int  RanV::bern(  double  p  ) 

{ 

double  u; 
u  =  mrandO  ; 

if  (  u  <=  p  )  return(l) ; 
else  return(O) ; 

}//end  function  bernO 

Function  binomO  :  Returns  binomial  R.V.  using  bernO. 

int  RanV::binom(  int  n,  double  p  ) 

{ 

int  sum=0 ; 

for  (  int  i=0;  i<n;  i++  ) 
sum  +=  bern(p) ; 

return  sum; 

}//end  function  binomO 
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Function  multnomO :  Returns  multinomial  R.  vector  using  binomQ .  Takes 
advantage  of  conditional  Xi  given  others  is  binomial. 

void  RanV : :multnom(  vector<double>  &mp,  vector<int>  &Xi  ) 

{ 

//initialize  variables  and  vectors 
int  m,  r,  j=l; 
m  =  (int)mp  [0] ; 
r  =  mp.sizeO  -  1; 

Xi . clear () ; 

for  (  int  i=0;  i<r;  i++  ) 

Xi . push_back(0) ; 

double  q=1.0; 

do{ 

Xi [ j  — 1]  =  binom(m,mp [j] /q) ; 
m  -=  Xi  [j-1]  ; 
q  -=  mp [ j ] ; 

j++; 

}while (  m  ! =  0  ) ; 

}//end  function  multnomO 

Function  ChiSquareO:  Returns  ChiSquare  R.V.  using  stdnormO  . 

double  RanV : : ChiSquare (  int  df  ) 

{ 

double  yLocal,  sum=0.0; 

for  (  int  i=0;  i<df;  i++  ){ 
yLocal  =  stdnormO  ; 
sum  +=  yLocal*yLocal ; 

}//end  for 

return  sum; 

}//end  function  ChiSquareO 

Function  gamma 0 :  Returns  gamma(a)  R.V.  using  GB  algorithm  (Cheng  1977) 
if  a  >  1 .  If  a  =  1  then  reduces  to  exponential.  If  a  <  1  then  uses 
rejection  algorithm  due  to  Ahrens  and  Dieter (1974) . 
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double  RanV::gamma(  double  alpha  ) 

{ 

if  (  alpha  <  1 . 0  ) 

return  gaml (alpha); 
else  if  (  alpha  ==  1.0  ) 

return  expn (alpha) ; 

else 

return  gam2 (alpha); 

}//end  function  gammaO 

//PRIVATE  FUNCTION 

double  RanV::gaml(  double  alphal  ) 

{ 

//See  Law  and  Kelton  for  description  of  algorithm 
double  b,  P,  Ul,  U2,  Y; 
double  Gaml=-1.0; 

b  =  (exp (1.0)  +  alphal)  /  exp (1.0); 
do  { 

Ul  =  mrandO  ; 

P  =  b  *  Ul; 

if  (  P  >  1.0  )  { 

Y  =  -log((b  -  P)/alphal); 

U2  =  mrandO  ; 

if  (  U2  <=  pow(Y,  alphal-1.0)  )  { 
Gaml  =  Y; 
break; 

}//end  if 

}else  { 

Y  =  pow(P,  1.0/alphal); 

U2  =  mrandO  ; 

if  (  U2  <=  exp(-Y)  ) 

Gaml  =  Y; 

}//end  if/else 
}while  (  Gaml  <  0.0  ); 

return  Gaml ; 

}//end  function  gaml() 

//PRIVATE  FUNCTION 

double  RanV::gam2(  double  alpha2  ) 

{ 

//See  Law  and  Kelton  for  description  of  algorithm 
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double  a,b,q,theta,d,ul ,u2; 
double  V , Y , Z , W ,  Gam2=-1.0; 

a  =  1.0/ (  sqrt (2 . 0*alpha2-l . 0)  ); 
b  =  alpha2  -  log(4.0); 
q  =  alpha2  +  1.0/a; 
theta  =  4.5; 
d  =  1.0  +  log(theta) ; 

if  (  alpha2  ==  0.0  ) 

Gam2  =  alpha2 ; 


do{ 

ul  =  mrandO  ; 
u2  =  mrandO  ; 

V  =  a*log(ul/ (1 . 0-ul) ) ; 

Y  =  alpha2*exp(V) ; 

Z  =  ul*ul*u2 ; 

W  =  b  +  q*V  -  Y; 

if  (  W  +  d  -  theta*Z  >=  0.0  ){ 

Gam2  =  Y; 
break; 

}//end  if 

if  (  W  >=  log(Z)  ){ 

Gam2  =  Y; 

}//end  if 

}while  (  Gam2  <  0.0  ); 

return  Gam2 ; 

}//end  function  gam2() 

Function  beta()  :  Returns  beta(al,a2)  R.V.  using  gamma () . 

double  RanV::beta(  double  al ,  double  a2  ) 

{ 

double  gl,  g2 ; 
gl  =  gamma(al) ; 
g2  =  gamma(a2) ; 

return (  gl/(gl+g2)  ); 

}//end  function  beta() 

Function  DirichletO:  Returns  Dirichlet  R. Vector  using  gammaO  . 
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void  RanV: : Dirichlet (  vector<double>  &Dp,  vector<double>  &XiD  ) 

{ 

int  k; 

k  =  Dp . size ()  ; 

vector<double>  yLocal; 

yLocal . clear () ; 

XiD.  clearO  ; 

for  (  int  i=0;  i<k;  i++  ){ 

XiD . push_back(0 . 0) ; 
yLocal. push_back(0. 0) ; 

}//end  for 

double  sum=0.0; 

for  (  int  i=0;  i<k;  i++  ){ 

yLocal [i]  =  gamma(Dp [i] ) ; 
sum  +=  yLocal [i] ; 

}//end  for 

for  (  int  i=0;  i<k;  i++  ) 

XiD[i]  =  yLocal [i] /sum; 

}//end  function  Dirichlet () 

Function  PoissonO :  Returns  Poisson  RV.  For  small  lamda  (<=12)  exploits 
relationship  with  exponential.  For  lamda  >  12,  uses  code  as  developed  in 
"Numerical  Recipes",  Press,  et  al . 

int  RanV: : Poisson (  double  lambda  ) 

{ 

const  double  PI  =  3.141592653589793238; 

static  double  sq,  alxm,  g,  oldl=(-1.0); 

//oldl  is  a  flag  for  whether  lamda  has  changed  since  last  call 

int  em; 

double  t,  tm,  y; 

if  (  lambda  <  12.0  )  {  //use  direct  method 

if  (  lambda  !=  oldl  )  { 
oldl  =  lambda; 

g  =  exp(  -lambda  );  //if  lambda  is  new,  compute  exponential 

}//end  if 
em  =  -1; 
t  =  1.0; 
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do  { 


++em; 

t  *=  mrandO  ; 

}while  (  t  >  g  ) ; 

}else  { 

if  (  lambda  !=  oldl  )  { 
oldl  =  lambda; 
sq  =  sqrt(  2.0*lambda  ); 
alxm  =  log(  lambda  ) ; 
g  =  lambda*alxm  -  gammln(lambda+l . 0) ; 

}//end  if 

do  -[ 

do  -[ 

y  =  tan(  PI*mrand()  )  ; 
tm  =  sq*y  +  lambda; 

}while  (  tm  <  0.0  )  ; 

em  =  (int)f loor (tm) ; 

t  =  0 . 9* (1 . 0+y*y) *exp(  em*alxm  -  gammln(em+l . 0)  -  g  ) 
}while  (  mrandO  >  t  )  ; 

}//end  else 

return  em; 

}//end  function  PoissonO 
//PRIVATE  FUNCTION 

//Returns  the  natural  log  of  the  Gamma  function  of  xx  for  xx>0 
double  RanV : : gammln (  const  double  xx  ) 

{ 

double  x,  y,  tmp,  ser; 

static  const  double  cof  [6]  =  {76.18009172947146,  -86.50532032941677, 
24.01409824083091,  -1.231739572450155,  0.001208650973866179, 
-0 . 000005395239384953} ; 

y  =  x  =  xx ; 
tmp  =  x  +  5.5; 
tmp  -=  (x+0 . 5) *log(tmp) ; 
ser  =  1.000000000190015; 
for  (  int  j=0;  j<6;  j++  ) 
ser  +=  cof[j]/++y; 

return  -tmp  +  log(  2 . 5066282746310005*ser/x  ); 

}//end  function  gammlnO 
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//DATA  HANDLING  FUNCTIONS 

Fuction  getdataQ :  Retrieves  data  from  file  and  stores  in  referenced 
vector.  Data  files  must  be  sequential  and  line,  space,  or  tab  delimited. 

void  RanV : : get data (  vector<double>  &data,  char*  filename  ) 

//open  input  file 

ifstream  inData(  filename,  ios::in  ); 

//exit  program  if  ifstream  could  not  open  file 
if  (  ! inData  )  { 

cerr  <<  "Data  file  could  not  be  opened"  <<  endl ; 
exit (  1  )  ; 

}//end  if 

//get  data 
data. clear () ; 
double  dvalue ; 

while  (  inData  >>  dvalue  ) 

data.push_back(  dvalue  ); 

//close  file 
inData. close () ; 

}//end  function  getdataO 

Fuction  cmomO :  Calculates  1st  and  2nd  central  moments/population 
mean  and  variance  of  vector  of  data. 

void  RanV::cmom(  vector<double>  &d,  double  &mean,  double  &var  ) 

{ 

//initialize  variables 
int  n  =  d.  sizeO  ; 
mean  =  0.0; 
var  =  0.0; 

//calculate  population  mean  and  variance 
for  (  int  i=0;  i<n;  i++  ) 
mean  +=  d  [i] ; 

mean  /=  n; 
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for  (  int  i=0;  i<n;  i++  ) 

var  +=  (d [i] -mean) * (d [i] -mean) ; 

var  /=  (n-1) ; 

}//end  function  cmomO 

Fuction  cov() :  Calculates  covariance  for  2  vectors  of  data. 

void  RanV::cov(  vector<double>  &dl,  vector<double>  &d2,  double  &cvar  ) 

{ 

//initialize  variables 

int  nl=dl . size () ,  n2=d2 . size () ,  n; 

n  =  min(nl ,n2) ; 

double  meanl=0.0,  mean2=0.0; 

//calculate  means  for  vectors 
for  (  int  i=0;  i<n;  i++  ){ 
meanl  +=  dl  [i]  ; 
mean2  +=  d2 [i] ; 

}//end  for 

meanl  /=  n; 
mean2  /=  n; 

//calculate  covariance 
for  (  int  i=0;  i<n;  i++  ) 

cvar  +=  (dl [i] -meanl) * (d2 [i] -mean2) ; 

cvar  /=  (n-1) ; 

}//end  function  cov() 

Fuction  quantO:  Calculates  specified  quantiles  for  vector  of  data. 
Supply  data  vector  and  vector  with  specified  quantiles  in  (0,1).  Returns 
sorted  data  vector;  quantile  values  in  place  of  passed  quantiles. 

void  RanV :: quant (  vector<double>  &dat ,  vector<double>  &q  ) 

//variables 
int  sd  =  dat.sizeQ; 
int  sq  =  q. size () ; 
int  loc; 
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//find  quantiles 

sort(  dat.beginO,  dat.endO  ); 

for  (  int  i=0;  i<sq;  i++  ){ 

loc  =  (int)ceil(  q[i]*sd  ); 
q[i]  =  dat[loc-l]; 

}//end  for 

}//end  function  quantO 
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Appendix  B 


Source  Code:  BayesRV.cpp 


The  following  is  verbatim  source  code  from  BayesRV.cpp  as  described 
in  chapters  2  and  3. 

//BayesRV.cpp,  Patrick  J.  Munson,  25  Apr  05 
//Member  function  definitions  for  class  BayesRV 
//generates  specialized  RVs  for  Bayesian  analysis 

#include  <iostream> 

#include  <cmath> 

#include  <vector> 

#include  <algorithm> 

#include  "RanV.h" 

using  namespace  std; 

//constructor 

BayesRV: : BayesRV (  int  seed  ) 

setSeed2(  seed  ); 

}//end  RanV  constructor 

//assigns  seed  value  to  stream2  for  use  in  RanV  objects 
void  BayesRV: : setSeed2(  int  s2  ) 

{ 

stream2  =  (  s2  >=  1  &&  s2  <=  10000  )  ?  s2  :  12; 

}//end  function  setSeed 


//PUBLIC  FUNCTIONS 

Function  TRexpnQ :  Returns  truncated  exponential (a)  random  variate  using 
truncated  uniform  and  the  inversion  method. 
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double  BayesRV: : TRexpn(  double  lambda,  double  a,  double  b  ) 

{ 

//variables  for  truncated  uniform 
double  TU,  Ua,  Ub; 

RanV  r(stream2) ;  //RanV  object  for  uniform  variates 

//Generate  Uniform  with  inverted  truncation  limits 

//We  can  remove  3  subtraction  operations  by  inverting  the  limits 

//and  subtracting  from  1 

Ua  =  exp(-lambda*b) ; 

Ub  =  exp(-lambda*a) ; 

TU  =  r . unif (Ua,Ub) ; 

//return  truncated  exponential 
return (  (-1.0/lambda)  *  log(TU)  ); 

}//end  function  TRexpnO 

Function  TRbetaQ :  Returns  truncated  Beta(a,b)  random  variate  using 
auxilliary  variable  technique  and  Gibbs  sampler. 

double  BayesRV: : TRbeta(  double  alpha,  double  beta,  double  a,  double  b  ) 

{ 

//variables 

double  X,  Y;  //X  given  Y;  Y  is  latent  variable  given  X 

double  xCL;  //truncation  limit  for  conditional  X 

//initialize  X 
X  =  alpha  /  (alpha+beta) ; 

RanV  r(stream2) ;  //RanV  object  for  uniform  variates 

//different  cases  for  beta=l,  beta<l,  and  beta>l 
if  (  beta  ==  1.0  )  { 

//simply  return  f(x)  =  a*x~(a-l) 

X  =  TBx(alpha,a,b) ; 

}else  if  (  beta  <  1.0  )  { 

//run  Gibbs  sampler  for  5  iterations 
for  (  int  i=0;  i<5;  i++  )  { 

//obtain  Y  given  X 

Y  =  r.unif(  0.0,(pow(  (1 . 0-X) , (beta-1 . 0)  ))  ); 

//now  obtain  X  given  Y 

xCL  =  1.0  -  pow(  Y, (1 . 0/ (beta-1 . 0) )  ); 
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xCL  =  max(  a,xCL  ); 

X  =  TBx (alpha, xCL, b) ; 

}//end  for 

}else  if  (  beta  >  1.0  )  { 

//run  Gibbs  sampler  for  5  iterations 
for  (  int  i=0;  i<5;  i++  )  { 

//obtain  Y  given  X 

Y  =  r.unif(  0.0,(pow(  (1 . 0-X) , (beta-1 . 0)  ))  ); 

//now  obtain  X  given  Y 

xCL  =  1.0  -  pow(  Y, (1 . 0/ (beta-1 . 0) )  ); 

xCL  =  min(  b,xCL  ); 

X  =  TBx (alpha, a, xCL) ; 

}//end  for 
}//end  if  else 

return  X; 

}//end  function  TRBetaO 

Function  TRgammaO :  Returns  truncated  gamma(a)  random  variate  using 
auxilliary  variable  technique  and  Gibbs  sampler. 

double  BayesRV: : TRgamma(  double  alpha,  double  a,  double  b  ) 

{ 

//variables 

double  X,  Y;  //X  given  Y;  Y  is  latent  variable  given  X 
double  bY ;  //to  store  new  truncation  limit  for  X  given  Y 

//initialize  X 
X  =  alpha; 

RanV  r(stream2) ;  //RanV  object  for  uniform  variates 

//run  Gibbs  sampler  for  5  iterations 
for  (  int  i=0;  i<5;  i++  )  -( 

//obtain  Y  given  X 
Y  =  r.unif(  0.0, (exp (  -X  ) )  ); 

//now  obtain  X  given  Y 
bY  =  —log (  Y  ); 
bY  =  min(  b,bY  )  ; 

//can  use  the  same  function  (TBx)  as  in  Beta 
X  =  TBx (alpha, a, bY) ; 

}//end  for 
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return  X; 

}//end  function  TRgamma 


Function  TRbiNormO :  Returns  truncated  bivariate  normal  random  vector 
using  auxilliary  variable  technique  and  Gibbs  sampler.  Pass  vector  for 
X  values,  mean  vector,  sigma  matrix,  and  truncation  matrix. 

void  BayesRV: :TRbiNorm(  double  X [] ,  double  Mu[] ,  double  Sig[] [2] ,  double  A [] [2]  ) 

{ 


//variables 
double  xlL [2] ,  x2L[2] ; 
double  B1 [2] ,  B2 [2] ; 
double  xMl ,  xM2; 
double  Sinv[2]  [2]  ; 
double  Sa,  Sb,  Sc,  Sd; 
double  Y,  Ya,  Yb; 

Yb  =  10000.0; 


//uniform  limits,  X  given  Y  is  U(a,b) 

//B  is  used  to  calculate  X  limits 
//stores  (X-mu)  to  simplify  calculations 
//for  storing  Sig  inverse 

//values  of  Sig  inverse — simplified  eqns 
//stores  Y  and  its  truncation  limits 

//right  limit  for  Y  given  X  is  infinity 


//initialize  X  and  xM 
X[0]  =  Mu  [0]  ; 

X  [1]  =  Mu  [1]  ; 
xMl  =  X[0]  -  Mu  [0]  ; 
xM2  =  X  [1]  -  Mu  [1]  ; 


//for  Damien(2001)  example  only 
A  [1]  [0]  =  0.5  -  sqrt  (  1- (X  [0] -0 . 5) * (X  [0] -0 . 5)  ); 

A  [1]  [1]  =  0.5  +  sqrt (  1- (X [0] -0 . 5) * (X [0] -0 . 5)  ); 

A  [0]  [0]  =  0.5  -  sqrt  (  1- (X  [1] -0 . 5) * (X  [1] -0 . 5)  ); 

A  [0]  [1]  =  0.5  +  sqrt  (  1- (X  [1] -0 . 5)  *  (X  [1] -0 . 5)  ); 


//calculate  the  inverse  of  the  Sig  matrix 
inv2x2(Sig,Sinv) ; 

Sa  =  Sinv [0] [0] ; 

Sb  =  Sinv [0] [1] ; 

Sc  =  Sinv  [1]  [0]  ; 

Sd  =  Sinv [1]  [1]  ; 

RanV  r(stream2);  //RanV  object  for  uniform  variates 

//run  Gibbs  sampler  for  5  iterations 
for  (  int  i=0;  i<5;  i++  )  { 
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//first  obtain  Y  given  X 

Ya  =  Sa*xMl*xMl  +  (Sb+Sc) *xMl*xM2  +  Sd*xM2*xM2; 

Y  =  TRexpn(0.5,Ya,Yb) ; 

//for  XI 

//solution  of  quadratic  eq  Y  >  (X-Mu) 1 [Siglnv] (X-Mu) 

B1 [0]  =  (1 . 0/ (2 . 0*Sa) )  *  (  -xM2* (Sb+Sc)  - 

sqrt (  4*Sa*(Y-xM2*xM2*Sd)  +  xM2*xM2* (Sb+Sc) * (Sb+Sc)  )  )  +  Mu[0] 
B1 [1]  =  (1 . 0/ (2 . 0*Sa) )  *  (  -xM2* (Sb+Sc)  + 

sqrt (  4*Sa*(Y-xM2*xM2*Sd)  +  xM2*xM2* (Sb+Sc) * (Sb+Sc)  )  )  +  Mu[0] 
//truncation  limits  for  X  given  Y 
xlL[0]  =  max(  A  [0]  [0]  ,B1  [0]  ); 
xlL  [1]  =  min(  A  [0]  [1]  ,B1  [1]  ); 

//now  obtain  XI 
X  [0]  =  r  .  unif  (xlL  [0]  ,xlL  [1] )  ; 

//recalculate  xMl 
xMl  =  X[0]  -  Mu  [0]  ; 

//for  Damien(2001)  example  only 
A  [1]  [0]  =  0.5  -  sqrt  (  1- (X  [0] -0 . 5) * (X  [0] -0 . 5)  ); 

A [1]  [1]  =  0.5  +  sqrt (  1- (X [0] -0 . 5) * (X  [0] -0 . 5)  ); 

//for  X2 

//solution  of  quadratic  eq  Y  >  (X-Mu) ’ Siglnv(X-Mu) 

B2 [0]  =  (1 . 0/ (2 . 0*Sa) )  *  (  -xMl* (Sb+Sc)  - 

sqrt (  4*Sa* (Y-xMl*xMl*Sd)  +  xMl*xMl* (Sb+Sc) * (Sb+Sc)  )  )  +  Mu[l] 
B2 [1]  =  (1 . 0/ (2 . 0*Sa) )  *  (  -xMl*(Sb+Sc)  + 

sqrt (  4*Sa* (Y-xMl*xMl*Sd)  +  xMl*xMl* (Sb+Sc) * (Sb+Sc)  )  )  +  Mu[l] 
//truncation  limits  for  X  given  Y 
x2L [0]  =  max(  A [1]  [0] ,B2  [0]  ); 
x2L  [1]  =  min(  A [1]  [1] ,B2[1]  ) ; 

//now  obtain  X2 
X  [  1]  =  r .  unif  (x2L  [0]  ,x2L  [1] )  ; 

//recalculate  xM2 
xM2  =  X  [1]  -  Mu  [1]  ; 

//for  Damien(2001)  example  only 
A  [0]  [0]  =  0.5  -  sqrt  (  1- (X  [1] -0 . 5) * (X  [1] -0 . 5)  ); 

A  [0]  [1]  =  0.5  +  sqrt  (  1- (X  [1] -0 . 5) * (X  [1] -0 . 5)  ); 


}//end  for 
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}//end  function  TRbiNorm 


//PRIVATE  FUNCTIONS 

Function  TBx():  Returns  X  given  Y  in  the  truncated  Beta  and  truncated 
Gamma  examples.  Here  F(x)  =  (x~alpha-a~alpha) / (b~alpha-a~alpha)  with 
X  on  the  set  (a,b) . 

double  BayesRV: :TBx(  double  &alpha,  double  &a,  double  &b  ) 

{ 

//variables 
double  U,  tmp; 

RanV  r(stream2);  //RanV  object 
U  =  r  .mrandQ  ; 

tmp  =  U* (pow(b , alpha) -pow (a, alpha) )  +  pow(a, alpha) ; 

//return  X,  given  latent  variable  Y,  for  truncated  Beta  and  Gamma 
return (  pow (tmp , (1 . 0/alpha) )  ); 

}//end  function  TBx() 

//calculate  and  return  the  inverse  of  a  2x2  matrix 
void  BayesRV: :inv2x2(  double  In[]  [2] ,  double  Out  []  [2]  ) 

{ 

double  detin;  //for  storing  the  determinant  of  In[] [] 

detin  =  ln[0]  [0]*ln[l]  [1]  -  ln[0]  [l]*In[l]  [0]  ; 

//calculate  the  inverse  and  return 
Out [0] [0]  =  In[l] [1] /detin; 

0ut[0][l]  =  -In  [0]  [1] /detin; 

0ut[l][0]  =  -In  [1]  [0] /detin; 

0ut[l][l]  =  ln[0]  [0] /detin; 

}//end  function  inv2x2() 
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Appendix  C 


Header  File:  RanV.h 


The  following  is  the  verbatim  header  file  RanV.h  as  described  in  chap¬ 


ter  2. 

//RanV.h 

//Includes  declarations  for  classes  RanV  and  BayesRV 
#ifndef  RANV_H 
#def ine  RANV_H 

//Declaration  of  class  RanV,  member  functions  are  defined  in  RanV.cpp 
class  RanV  { 

public : 


//DEFAULT  CONSTRUCTOR 
RanV(  int  =  12  )  ; 

//set  seed 

void  setSeedC  int  ) ; 

//UNIFORM 
double  mrandO  ; 

double  unif (  double  a,  double  b  ); 

void  mrandst (  std: : vector<double>  &,  int  str  ); 

void  mrandgt (  std: : vector<double>  &,  int  str  ); 

//NONUNIFORM 
double  stdnormO  ; 

double  norm(  double  mu,  double  sig2  ) ; 

double  logn(  double  mu,  double  sig2  ) ; 

double  expn(  double  lambda  ) ; 

double  weibull(  double  shape,  double  scale  ); 

int  bern(  double  p  ) ; 

int  binom(  int  n,  double  p  ) ; 

void  multnom(  std: : vector<double>  &mp,  std: : vector<int>  &Xi  ); 
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double  ChiSquare(  int  df  ); 
double  gamma (  double  alpha  ) ; 
double  beta(  double  al,  double  a2  ); 

void  Dirichlet(  std: : vector<double>  &Dp,  std: : vector<double>  &XiD  ); 
int  Poisson (  double  lambda  ) ; 

//DATA  HANDLING 

void  getdata(  std: : vector<double>  &,  char*  ); 

void  cmom(  std : : vector<double>  &,  double  &mean,  double  &var  ); 

void  cov(  std : : vector<double>  &,  std : : vector<double>  &,  double  &cvar  ); 

void  quant (  std: : vector<double>  &,  std : : vector<double>  &  ); 


private : 


//Seed  variable 
int  stream; 

//Gamma  functions:  depends  on  alpha 
double  garni (  double  alphal  ) ; 
double  gam2(  double  alpha2  ); 

//for  use  in  PoissonO 

double  gammln(  const  double  xx  ) ; 


};//end  class  RanV 

//Declaration  of  class  BayesRV,  member  functions  are  defined  in  BayesRV.cpp 
class  BayesRV  { 

public : 

BayesRV (  int  =  12  )  ; 
void  setSeed2(  int  ); 


double  TRexpn(  double  lambda,  double  a,  double  b  ); 

double  TRbeta(  double  alpha,  double  beta,  double  a,  double  b  ) ; 

double  TRgamma(  double  alpha,  double  a,  double  b  ); 

void  TRbiNorm(  double  X [] ,  double  Mu[] ,  double  Sig  []  [2] ,  double  A []  [2]  ); 


private : 


int  stream2; 

double  TBx(  double  &,  double  &,  double  &  ); 
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void  inv2x2(  double  []  [2] ,  double  []  [2]  ); 


};//end  class  BayesRV 
#endif 
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